      SUBROUTINE QMC(g0,g,xmu,dtau,u,iscf,idum,acc,nrat,nswpIN,nn)         
! ----------  ----------  ----------  ----------  ----------  ----------
!     IMPLICIT REAL*8 (A-H,O-Z)
      include 'param.dat'

! ----------  ----------  ----------  ----------  ----------  ----------
#ifdef MPI
      include 'mpif.h'
#endif
! ----------  ----------  ----------  ----------  ----------  ----------

      REAl*8                                                      
     &   g0(L,Nd),   g(L,Nd)                                                
     & , ga(L,L,Nd), gx(L,L,Nd), gst(L,L,Nd)                                
     & , v(L),       vn(L,Nf)                                               
     & , del(L,L)                                                           
     & , xg(-L:L),   xgf(-L:L),  G00(-L:L),  xs(0:L)                        
     & , xf(Nd,Nf)                                                          
     & , u(Nf),      xlam(Nf)                                               
     & , DOcc(Nf),   nn(Nf)                                                 
      REAl*8  xmud
!                                                                           
      DIMENSION s(2*L,1,Nf) ! it is to save memory space                    
      INTEGER ipoint(Nf,2)                                                  
      INTEGER myid,ierr,numprocs,master                                     
!     COMMON/IFs/   IF1,IFf                                                 
!                                                                           
      COMMON/varmpi/myid,ierr,numprocs,master                               
      COMMON/ising/ vn                                                      
                                                                            
! ----------  ----------  ----------  ----------  ----------  ----------    
	IF1=1                                                               
	IFF=1                                                               
                                                                            
!      PRINT *, 'IF1, IFF=',  IF1, IFF                                      
!      PRINT *, 'myid,ierr,numprocs,master', myid,ierr,numprocs,master      
	 xmud=dtau*xmu                                                      
                                                                            
      IF(idum.eq.0) idum=123457                                             
      nsweep=nswpIN/numprocs                                                
      polar=0.5d0                                                           
      ndirty=100                                                            
      ncor=2                                                                
      IF(iscf.eq.1) THEN                                                    
        nwarm=100                                                           
      ELSE                                                                  
        nwarm=30                                                            
      ENDIF                                                                 
      DO in=1,Nf                                                            
    	  DOcc(in)=0.d0                                                     
  	    nn(in)=0.d0                                                     
      ENDDO                                                                 
! ----------  ----------  ----------  ----------  ----------  ----------    
      DO j=1,L                                                              
         DO  i=1,L                                                          
            del(i,j)=0.d0                                                   
            del(i,i)=1.d0                                                   
         ENDDO                                                              
      ENDDO                                                                 
!                                                                           
      DO kc=1,Nd                                                            
       DO j=1,L                                                             
         DO i=1,L                                                           
            gst(i,j,kc)=0.d0                                                
         ENDDO                                                              
        ENDDO                                                               
      ENDDO                                                                 
!                                                                           
! initialize: xf                                                            
! xf contains in the sign of the interaction (ie. +=ferro                   
! and -=antIFerro) given the channel (1-Nd) and the pseuDO spin (1-Nf)      
	DO i=1,Nd                                                           
	 DO j=1,Nf                                                          
	  xf(i,j)=0.d0                                                      
	 ENDDO                                                              
	ENDDO                                                               
!                                                                           
!       Sigma(i,ij)  =  1 IF i<j                                            
!                    = -1 IF i>j                                            
        ij=0                                                                
        DO i=1,Nd-1                                                         
         DO j=i+1,Nd                                                        
          ij=ij+1                                                           
          xf(i,ij)=1.d0                                                     
          xf(j,ij)=-1.d0                                                    
         ENDDO                                                              
        ENDDO                                                               
!                                                                           
! initialization: ipoint(:2)                                                
! first entry (1-Nf) is which pseuDO spin flips                             
! and the second entry (1-2) corresponds to                                 
! 1 to ferro-coupling (+), and 2 to antIFerro (-)                           
        ij=0                                                                
        DO i=1,Nd-1                                                         
         DO j=i+1,Nd                                                        
          ij=ij+1                                                           
          ipoint(ij,1)=i                                                    
          ipoint(ij,2)=j                                                    
         ENDDO                                                              
        ENDDO                                                               
!                                                                           
	DO i=1,Nf                                                           
          z=dtau*u(i)/2.d0                                                  
          z=exp(z)                                                          
          xlam(i)=dlog(z+dsqrt(z**2-1.d0))                                  
	ENDDO                                                               
!     initialization: ranDOm numbers generator                              
      r=ranw(idum)                                                          
!// Minus - change for QMC G-notation //                                    
	DO kc=1,Nd                                                          
         DO i=0,L-1                                                         
          g00(i)=-g0(i+1,kc)                                                
         ENDDO                                                              
         DO i=1,L-1                                                         
          g00(-i)=-g00(L-i)                                                 
         ENDDO                                                              
!?      g00(-L)=-g00(0)                                                     
	 DO  ir=1,L                                                         
	  DO  is=1,L                                                        
           gx(is,ir,kc)=g00(is-ir)                                          
	  ENDDO                                                             
	 ENDDO                                                              
	ENDDO                                                               
!                                                                           
      IF(iscf.eq.1) THEN                                                    
	DO iv=1,Nf                                                          
	 xlam0=xlam(iv)                                                     
	 CALL initial(xlam0,v,polar,idum)                                   
	 DO i=1,L                                                           
	  vn(i,iv)=v(i)                                                     
	 ENDDO                                                              
	ENDDO                                                               
      ENDIF                                                                 
!                                                                           
!     IF( myid.eq.master) THEN                                              
!        print*,'Ising fields - start'                                      
!       DO iv=1,Nf                                                          
!        write(*,'(8f10.3)')(vn(i,iv),i=1,L)                                
!       ENDDO                                                               
!       DO kc=1,Nd                                                          
!       print*,'G(i)=',kc                                                   
!       write(6,'(8f10.4)')((gx(i,j,kc),i=1,L),j=1,L)                       
!       ENDDO                                                               
!     ENDIF  ! master                                                       
!                                                                           
	DO kc=1,Nd                                                          
          CALL gnewclean(ga,xmud,vn,gx,del,xf,kc)                           
	ENDDO                                                               
! ----------  ----------  ----------  ----------  ----------  ----------    
! Main loop DOes nsweeps: a clean update comes after ndirty dirty updates   
!                                                                           
! kx= # of measurements                                                     
! irr= # of accepted flips (after the warm-up)                              
! nrat= # of negative det encountered                                       
! ----------  ----------  ----------  ----------  ----------  ----------    
	kxmax=0                                                             
	kx=0                                                                
      	irr=0                                                               
      	nrat=0                                                              
!                                                                           
      	DO 2 k=1,nsweep                                                     
         kk=mod(k,ndirty)                                                   
         kcor=mod(k,ncor)                                                   
         DO 5 j=1,L                                                         
	  DO 55 iv=1,Nf                                                     
            dv=2.d0*vn(j,iv)                                                
!// calculates the determinant ratio //                                     
            kup=ipoint(iv,1)                                                
	    kdw=ipoint(iv,2)                                                
            ratup=1.d0+(1.d0-ga(j,j,kup))*(dexp(-dv)-1.d0)                  
            ratdw=1.d0+(1.d0-ga(j,j,kdw))*(dexp(dv)-1.d0)                   
            rat=ratup*ratdw                                                 
!                                                                           
            IF(rat.lt.0.)THEN                                               
               nrat=nrat+1                                                  
            END IF                                                          
            rat=rat/(1.d0+rat)                                              
!// for num rec                                                             
	    r=ranw(idum)                                                    
            IF(rat.gt.r)THEN                                                
               IF(k.gt.nwarm)THEN                                           
                  irr=irr+1                                                 
               END IF                                                       
               vn(j,iv)=-vn(j,iv)                                           
               IF(kk.eq.0)THEN                                              
      		  CALL gnewclean(ga,xmud,vn,gx,del,xf,kup)                  
      		  CALL gnewclean(ga,xmud,vn,gx,del,xf,kdw)                  
               ELSE                                                         
                  CALL gnew(ga,vn,j,iv,del,1.d0,kup)                        
                  CALL gnew(ga,vn,j,iv,del,-1.d0,kdw)                       
               ENDIF                                                        
            ENDIF                                                           
55  	  CONTINUE                                                          
5       CONTINUE                                                            
!// store the measurements //                                               
	IF(kcor.eq.0.and.k.gt.nwarm)THEN                                    
	kx=kx+1                                                             
	DO kc=1,Nd                                                          
!// store GF //                                                             
	  DO jx=1,L                                                         
            DO ix=1,L                                                       
              gst(ix,jx,kc)=gst(ix,jx,kc)+ga(ix,jx,kc)                      
     	    ENDDO                                                           
	  ENDDO                                                             
	ENDDO                                                               
!                                                                           
! ----------  ----------  ----------  ----------  ----------  ----------    
	  ij=0                                                              
        DO ik=1,2*Nlm-1                                                     
         DO jk=ik+1,2*Nlm                                                   
		ij=ij+1                                                     
          DO ix=1,L                                                         
!                                                                           
           gak  =1.-ga(ix,ix,ik)                                            
           gakN =1.-ga(ix,ix,jk)                                            
!                                                                           
           DOcc(ij)=DOcc(ij)+gak*gakN                                       
!                                                                           
          ENDDO                                                             
         ENDDO                                                              
        ENDDO                                                               
! ----------  ----------  ----------  ----------  ----------  ----------    
	IF(IF1.NE.1) THEN                                                   
!// store the Ising spins //                                                
        DO iq=1,L                                                           
 	  DO iv=1,Nf                                                        
          s(iq,kx,iv)=vn(iq,iv)/abs(vn(iq,iv))                              
          s(iq+L,kx,iv)=s(iq,kx,iv)                                         
	  ENDDO                                                             
        ENDDO                                                               
 	ENDIF ! (IF1.NE.1)                                                  
!                                                                           
 	ENDIF                                                               
!                                                                           
 2    CONTINUE  ! END of the main QMC loop                                  
! ----------  ----------  ----------  ----------  ----------  ----------    
!     IF( myid.eq.master) THEN                                              
!        print*,'Ising fields - END QMC'                                    
!       DO iv=1,Nf                                                          
!        write(*,'(16f5.0)')(vn(i,iv),i=1,L)                                
!       ENDDO                                                               
!     ENDIF  ! master                                                       
!// acceptance rate :                                                       
      acc=dfloat(irr)/dfloat((nsweep-nwarm)*L*Nf)                           
!// normalization of GF //                                                  
!      kxmax=kx*numprocs                                                    
      kxmax=kx                                                              
      DO kc=1,Nd                                                            
	DO jx=1,L                                                           
          DO ix=1,L                                                         
            gst(ix,jx,kc)=gst(ix,jx,kc)/DFLOAT(kxmax)                       
	  ENDDO                                                             
	ENDDO                                                               
      ENDDO                                                                 
!                                                                           
!// and the DOuble occupation //                                            
!     IF(myid .eq. master) THEN  ! master                                   
!       print*,'i, u(i),DOuble occ(i)='                                     
!       IF(IF.eq.0)THEN                                                     
!                                                                           
         DO k=1,Nf                                                          
            nn(k)=DOcc(k)/DFLOAT(kxmax*L*numprocs)                          
         ENDDO                                                              
!                                                                           
!       ENDIF                                                               
!     ENDIF !master                                                         
! ----------  ----------  ----------  ----------  ----------  ----------    
!// wrap-around //                                                          
	DO kc=1,Nd                                                          
!// negative times (-) //                                                   
          DO j=0,L-1                                                        
            xg(-j)=0.                                                       
            DO i=1,L-j                                                      
              xg(-j)=xg(-j)+gst(i,i+j,kc)                                   
            ENDDO                                                           
          ENDDO                                                             
!// positive times (+) //                                                   
          DO j=1,L-1                                                        
            xg(j)=0.                                                        
            DO i=1,L-j                                                      
              xg(j)=xg(j)+gst(i+j,i,kc)                                     
            ENDDO                                                           
          ENDDO                                                             
!// it is enforced the time "antiperiodicity" //                            
        DO i=1,L-1                                                          
          xgf(i)=(xg(i)-xg(i-L))/dfloat(L)                                  
          xgf(i-L)=-xgf(i)                                                  
        ENDDO                                                               
        xgf(0)=xg(0)/dfloat(L)                                              
        xgf(-L)=-xgf(0)                                                     
!// change MINUS sign from QMC //                                           
        DO i=0,L-1                                                          
          g(i+1,kc)=-xgf(i)                                                 
	ENDDO                                                               
!                                                                           
	ENDDO  ! kc=1,Nd                                                    
! ----------  ----------  ----------  ----------  ----------  ----------    
!// computation of susceptibilities //                                      
      IF(IFf.NE.1) THEN                                                     
      IF(IF1.eq.0)THEN                                                      
!// get chis: first from the ising fields //                                
!                                                                           
      DO iv=1,Nf                                                            
!!!!!have to fix this!!!                                                    
! one has to have Nf dIFferent suceptibilities!!!!                          
        kxa=1                                                               
        xxx=1.                                                              
        xsa=0.                                                              
        xs2=0.                                                              
        DO  ii=0,L                                                          
          xs(ii)=0.                                                         
        ENDDO                                                               
        DO iq=1,kxa                                                         
          DO itq=0,L-1                                                      
            DO il=1,L                                                       
              itk=il+itq                                                    
              xs(itq)=xs(itq)+s(itk,iq,iv)*s(il,iq,iv)                      
              xsa=xsa+s(il,iq,iv)                                           
              xs2=xs2+s(il,iq,iv)**2                                        
            ENDDO                                                           
          ENDDO                                                             
        ENDDO                                                               
        DO ii=1,L-1                                                         
          xs(ii)=xxx*xs(ii)/dfloat(kxa*L)                                   
          write(43,*)dfloat(ii)/dtau/dfloat(L),xs(ii)                       
        ENDDO                                                               
      ENdDO                                                                 
!                                                                           
      ENDIF ! (IF1.eq.0)                                                    
      ENDIF ! (IFf.NE.0)                                                    
! ----------  ----------  ----------  ----------  ----------  ----------    
      RETURN                                                                
      END                                                                   
! ----------  ----------  ----------  ----------  ----------  ----------    
! ----------  ----------  ----------  ----------  ----------  ----------    
!// dirty update of GF //                                                   
!                                                                           
      SUBROUTINE gnew(g,vn,j,iv,del,xflag,kc)                               
      include 'param.dat'                                                   
      REAl*8 g(L,L,Nd),vn(L,Nf),del(L,L),d(L,L)                   
!                                                                           
        dv=xflag*2.*vn(j,iv)                                                
        ee=exp(dv)-1.d0                                                     
        a=ee/(1.d0+(1.d0-g(j,j,kc))*ee)                                     
        DO i2=1,L                                                           
          DO i1=1,L                                                         
            d(i1,i2)=g(i1,i2,kc)+(g(i1,j,kc)-del(i1,j))*(a*g(j,i2,kc))      
          ENDDO                                                             
        ENDDO                                                               
        DO i2=1,L                                                           
          DO i1=1,L                                                         
            g(i1,i2,kc)=d(i1,i2)                                            
          ENDDO                                                             
        ENDDO                                                               
      RETURN                                                                
      END                                                                   
! ----------  ----------  ----------  ----------  ----------  ----------    
! ----------  ----------  ----------  ----------  ----------  ----------    
!// clean update of GF //                                                   
!                                                                           
      SUBROUTINE gnewclean(g,xmud,vn,gx,del,xf,kc)                          
      include 'param.dat'                                                   
      REAl*8  g(L,L,Nd), vn(L,Nf)                                 
     & , b(L,L), binv(L,L), del(L,L)                                        
     & , gx(L,L,Nd), ee(L)                                                  
     & , xf(Nd,Nf)                                                          
!                                                                           
	DO  i=1,L                                                           
	vv=0.d0                                                             
	 DO iv=1,Nf                                                         
	  vv=vv+xf(kc,iv)*vn(i,iv)                                          
	 ENDDO                                                              
	  vv=vv+xmud                                                        
 	ee(i)=dexp(vv)-1.d0                                                 
	ENDDO                                                               
!                                                                           
      DO j=1,L                                                              
         DO  i=1,L                                                          
            b(i,j)=del(i,j)-ee(j)*(gx(i,j,kc)-del(i,j))                     
	 ENDDO                                                              
      ENDDO                                                                 
!                                                                           
      CALL inverse(b,binv)                                                  
!                                                                           
      DO i1=1,L                                                             
         DO i2=1,L                                                          
            xdum=0.d0                                                       
            DO i=1,L                                                        
               xdum=xdum+binv(i1,i)*gx(i,i2,kc)                             
	    ENDDO                                                           
           g(i1,i2,kc)=xdum                                                 
          ENDDO                                                             
       ENDDO                                                                
!                                                                           
      RETURN                                                                
      END                                                                   
! ----------  ----------  ----------  ----------  ----------  ----------    
! ----------  ----------  ----------  ----------  ----------  ----------    
!// initialize the vector v of Ising fields //                              
!                                                                           
      SUBROUTINE initial(xlam,v,polar,idum)                                 
      include 'param.dat'                                                   
      REAl*8 v(L)                                                 
!                                                                           
      DO i=1,L                                                              
         s=1.d0                                                             
	  r=ranw(idum)                                                      
         IF(r.gt.polar) s=-1.d0                                             
         v(i)=xlam*s                                                        
      ENDDO                                                                 
      RETURN                                                                
      END                                                                   
! ----------  ----------  ----------  ----------  ----------  ----------    
! ----------  ----------  ----------  ----------  ----------  ----------    
